Investigating the migration of immiscible contaminant fluid flow in homogeneous and heterogeneous aquifers with high-precision numerical simulations

Numerical modeling of the migration of three-phase immiscible fluid flow in variably saturated zones is challenging due to the different behavior of the system between unsaturated and saturated zones. This behavior results in the use of different numerical methods for the numerical simulation of the fluid flow depending on whether it is in the unsaturated or saturated zones. This paper shows that using a high-resolution shock-capturing conservative method to resolve the nonlinear governing coupled partial differential equations of a three-phase immiscible fluid flow allows the numerical simulation of the system through both zones providing a unitary vision (and resolution) of the migration of an immiscible contaminant problem within a porous medium. In particular, using different initial scenarios (including impermeable “lenses” in heterogeneous aquifers), three-dimensional numerical simulation results are presented on the temporal evolution of the contaminant migration following the saturation profiles of the three-phases fluids flow in variably saturated zones. It is considered either light nonaqueous phase liquid with a density less than the water, or dense nonaqueous phase liquid, which has densities greater than the water initially released in unsaturated dry soil. Our study shows that the fate of the migration of immiscible contaminants in variably saturated zones can be accurately described, using a unique mathematical conservative model, with different evolution depending on the value of the system’s physical parameters, including the contaminant density, and accurately tracking the evolution of the sharp (shock) contaminant front.


Introduction
Multiphase flow problems refer to the simultaneous flow of two or more fluids separated by sharp interphases. They are observed in natural phenomena like multiphase flow within a porous structure, blood flow, nuclear reactions, oil and gas industry applications, etc. For these reasons, suitable numerical models are necessary to predict their physical behavior accurately. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 However, modeling multiphase fluid flow is challenging since a universally robust and accurate solution methodology has not yet been identified. It is important to mention that this work is only concerned with unsteady three-phase fluids flow in the variably saturated zone.
A multiphase fluid flow dynamics are governed by coupled conserved partial differential equations for each fluid flow, based on the Darcy Law and the mass and momentum conservation principles. They are written as a function of saturation, capillary pressure, permeability, porosity, density, and viscosity of each fluid flow. Since the capillary pressure and permeability of each phase is a function of the saturation, these equations are highly nonlinear and, due to the gravity and pressure gradients, are responsible for creating sharp (shocks) front and rarefaction, which may introduce significant errors in the numerical simulations. For this reason, it is not easy to obtain a numerical solution that converges to the physical solution, although there is an increasing effort to construct conservative numerical solution methods. See, for example, recent reviews for the numerical resolution of the Richards' equation for variably saturated flow [1,2], immiscible fluid flow in unsaturated zones [3], transport modeling in the heterogeneous porous medium [4,5], including nonaqueous phase liquid (NAPL) [6,7], and variably saturated zones [8][9][10]. Several two-phase flow models exist in the literature that uses high-resolution numerical schemes, such as high-resolution central upwind scheme [11]. They have been used as a numerical approximation to solve the two-phase fluid flow model [12]. Or the numerical solution of hyperbolic conservation laws such as a space-time Conservation Element-Solution Element (CESE) high-resolution scheme for computing the transport of a diffusing pollutant in shallow flows [13], using a Kinetic Flux Vector Splitting scheme [14], using a fifth-order Weighted Essentially Non-Oscillatory (WENO) scheme [15]. Finite volume WENO scheme and discontinuous Galerkin method have been used to solve multiphase flow models [16][17][18][19][20], a second-order accurate difference method for nonlinear conservation laws in three-dimensions [21], a shock-capturing numerical scheme for compressible multicomponent problems [22] and multi-medium flows [23]. Other examples of two-phase fluid flow applications for oil-water fluids flow can be found in [24,25].
The method used in this paper has been introduced in Ref. [26]. It is based on the high-resolution shock-capturing flux (HRSC) conservative method [11,27,28] to follow sharp discontinuities accurately and temporal dynamics of three-phase immiscible fluid flow in a porous medium. Several validation tests were performed in [26] to verify the accuracy of the HRSC method and the CactusHydro code. These tests include a comparison with an analytical model such as the Burgers' equation, the Buckley-Leverett model and a comparison with a twodimensional unsaturated-saturated water flow model [10], together with the sand tank experimental data [29]. They show the absence of spurious oscillations in the solution and convergence to the "weak" solution as the grid is refined. The time evolution is performed using a forward in time explicit method rather than the most used, implicit one in which the discretization is based on a "backward in time" evolution. That requires the time step to be sufficiently small since the method is "conditionally stable". The implicit methods, in contrast, are "unconditionally stable" but very expensive from the computational point of view and may lead to mass balance errors. See for example [30,31] for two-phase flow using an implicit pressure and explicit saturation schemes, three-phase fluid flow [32], finite element methods [33][34][35][36].
This work shows numerical results on the fate of either light NAPL or dense NAPL initially released on a dry unsaturated zone that migrates, depending on its density to the saturated aquifer one. In the literature, these scenarios are treated very differently. Here different initial scenarios are investigated (including impermeable "lenses" in heterogeneous aquifers). It is shown that using a unique mathematical model that uses an HRSC method, it is possible to numerical simulate the system in both variable saturated zones, giving a unitary vision and numerical resolution of the migration of immiscible contaminants in a porous medium. In particular, it follows the temporal (and three-dimensional spatial) evolution of the saturation profile of the three-phases fluid along the variably saturated zone and their conservation law all over the simulation.

Governing equations of three-phase immiscible fluid flow
This section briefly describes the numerical model introduced in Ref. [26]. Consider a threedimensional three-phase fluid flow in a porous medium composed of nonaqueous (n), water (w), air (a), and a variably saturated zone. Using the conservation equation for the mass and momentum for each fluid phase together with the Darcy's velocity for each phase, we obtain the governing coupled partial differential equations (PDEs) for each phase-fluid flow, where x i = (x, y, z) are the spatial cartesian coordinates, and t is the time coordinate, ρ α is the density M , k rα is the dimensionless relative permeability of phase α, k ij is the absolute permeability tensor [L 2 ], g denotes the gravitational acceleration L , ϕ is the porosity, S α is the dimensionless volumetric saturation of phase α which satisfy the relation, In Eqs (1)-(3) it is used the pressure p a (that is the air pressure when S a is different from zero), and it is considered the capillary pressure for the air-water phase, p caw = (p a − p w ), and the capillary pressure for the air-nonaqueous phase, p can = (p a − p n ), where it is substituted p w = p a − p caw , and p n = p a − p can . The third capillary pressure, the nonaqueous-water phase can be deduced from the other two, and is given by, p cnw = (p n − p w ) = (p caw − p can ) (in contrast to Refs. [37,38], where the air gradient pressure is assumed negligible). The relative permeabilities k rw , k rn and k ra and the capillary pressures are function of saturations, k rα = k rα (S a , S n , S w ), p can = p can (S a , S n , S w ) and p caw = p caw (S a , S n , S w ). This paper uses the van Genuchten model [39,40] for the air-water and the air-nonaqueous capillary pressure (see subsection 2.3). But it is worth noticing that different solutions choices correspond to different porous mediums. The numerical solution and the method applied here are not affected by any particular choice. The absolute permeability k ij depends on the properties of the porous medium. The porosity ϕ is a function of pressure and can be linearly approximated to, ϕ = ϕ 0 [1 + c R (p − p 0 )], where c R is the rock compressibility, ϕ 0 is the porosity at p 0 , which is considered the atmospheric, and p is the pressure (that will be associated with p a ).
It is convenient to define the product of the porosity ϕ and the saturation for each phase as, σ w � ϕS w , σ n � ϕS n , σ a � ϕS a . Then, Eq (4) can be written as, The left-hand side of (1)-(3), assuming constant density-viscosity for each phase, becomes @s a @t þ @s n @t þ @s w @t where @s a do not depend on the spatial derivative of the saturation, while depends on the spatial derivative of the saturation. The PDEs system to be numerically resolved is composed by Eqs (5) and (6) and the variables p, σ w , σ n , σ a together with the functional form of the relative permeabilities and capillary pressures, written in terms of the saturations (see subsection 2.3).

High-resolution shock-capturing conservative numerical method
The system of PDEs (5) and (6) governing a three-phase immiscible fluid is naturally formulated in terms of conservation laws. In this system, the dominant part (in the regime of interest) is the hyperbolic one (Eq (7)) and is responsible for forming shocks and propagating the density fronts. Indeed, to accurately reproduce the dynamics of the density discontinuity in time, one needs to separate the methods used to discretize the propagation (hyperbolic) part from the diffusive (parabolic) part (Eq (8)). This separation is of fundamental importance since explicit methods in time evolution should be preferred when the propagation dominates over the diffusion. In contrast, implicit (in time) methods should be preferred when diffusion dominates over propagation. The importance of using conservative formulations, i.e., methods based on conservation law, is due to two fundamental theorems. The first one is due to Lax and Wendorff [27] and the second one by Hou and LeFloch [28]. Namely, conservative numerical schemes, if convergent, do converge to the weak solution of the problem. The second theorem states that non-conservative numerical schemes do not converge to the correct solution if a shock wave (or discontinuity) is present in the flow. These two theorems state that if a conservative formulation is used, then it is guaranteed that the numerical solution will converge to the correct solution. On the contrary, if a conservative formulation is not used, the numerical solution will converge to the incorrect solution when the flow develops a discontinuity.
There is a vast literature on numerical methods that take advantage of the conservative nature of the equation and are referred to as HRSC methods that accurately reproduce the discontinuous features of the solutions. For the system studied in this paper, information on the characteristic is not explicitly available. Indeed, we are using a variant of the HRSC central schemed discussed in Kurganov and Tadmor (KT) [11] that also has the advantage of dealing with a diffusive part in the conservative equations that is usually assumed not present in other schemes. A comprehensive and almost complete discussion of most of the possible HRSC methods can be found in Ref. [41] dedicated, in a different context, to the study of Relativistic Hydrodynamics. Here is briefly described the details of the method in one spatial dimension and one component fluid. Its extension to more components and three dimensions is straightforward. The equation to be solved is, where F(u(x)) is the flux associated with the hyperbolic part, and Q(u(x), du(x)/dx) is the parabolic part of the flux. Notice that it is not explicitly written the dependence on t in all quantities. How this decomposition is achieved for the system in discussion is shown in (7,8). The variable is discretized assuming, at any time, a given value (u i ) on each grid point (x i ). For intermediate points is assumed a piecewise-linear reconstruction inside each cell, where s i is the "slope" of the linear reconstruction, and they are constructed using total-variation-diminishing (TVD) properties (needed for the theorem named above) and based on the min-mod slope limiter [Kolgan 72, val Leer 79]. Indeed, the variable is assumed to have a jump discontinuity at the cell border, i.e., different values on the right and the left at the points x i+1/2 called u þ i and u À iþ1 . This implies that the computed values of the flux F(u(x)) and Q(u(x), du(x)/dx) = Q(u(x), u 0 (x)) will be different at the cell boundary but these values may be used to compute the effective flux at the cell boundary H j+1/2 and P j+1/2 . Here H j+1/2 is . Once a prescription to construct this flux is given, one has a conservative flux method for the solution, namely one has just to numerical solve in time the ordinary differential equations (Methods of Lines) associated to each grid point: The KT method amount to the following choice for the numerical fluxes: where the a j+1/2 depends on the explicit functional dependence of the flux F(x, u(x)) that, in this case, is not explicitly known without assuming an explicit dependence on the permeabilities and their derivatives. However, independently on the choice of the permeabilities, the flux on the two sides always have the same sign and indeed, to use for the H j+1/2 the values computed on the "left" or on the "right" depend on the sign of F(x, u(x)) at the interface between two adjacent cells, namely The Kurganov-Tadmor (KT) scheme avoids the local Riemann problem. The method belongs to the Monotonic Upstream-centered Scheme for Conservation Law (MUSCL) class suggested by van Leer in 1973. The KT scheme archives second-order accuracy in space using analytical information on the flux form. However, the variation of the method used here that uses the first-order upwind formula for the fluxes and the min-mod flux limiter avoids using any analytical information (besides monotonicity) and, indeed, only the point values of the flux are required. This is of great advantage since it can use tabulated values for permeabilities. The only penalty is that it is a first-order accuracy at the discontinuities and over the whole simulation grid. This is not a big drawback since one must keep in mind that any method will be first-order accuracy at the physical discontinuities, and the method is targeted to study and follow the evolution of sharp discontinuity fronts.

Three permeabilities and capillary pressures model
The relative permeabilities for three phases are given in Ref. [39] and listed below, where S et is the total effective liquid saturation defined in terms of the irreducible wetting phase saturation S wir . This paper uses the van Genuchten model [38] where the effective saturation, S e , is given by, ð Þ , and α and n are model parameters. It can be solved for p c , where n ¼ 1 1À m , and p c0 = α −1 is the capillary pressure at S e = 0. Since p caw = p can + p cnw , the capillary pressures are given by,

Three-dimensional immiscible fluid flow numerical simulations results
This section presents the results of seven different three-dimensional numerical simulations "real" examples of an immiscible contaminant released to the environment in the unsaturated zone that migrates toward the saturated (aquifer) under the effect of the gravity force; exploring different scenarios in both homogeneous and heterogeneous aquifer systems. It is considered a three-phase fluid model composed of water, air, and nonaqueous phase liquid (light or dense NAPL) and investigates the temporal evolution of the migration of the contaminant following the saturation profiles of each phase along the variably saturated zone. We assumed constant density and constant viscosity for each phase. The effects of volatilization/biodegradation and dissolution are not considered. The dimension of the grid has been chosen to cover the minimum possible part of the aquifer but big enough to avoid the boundaries effect (finite grid size) on the dynamics of the contaminant. The porous medium is composed of an unsaturated dry zone (air-NAPL) and a saturated one (filled up with water, colored in blue) separated by the groundwater table located at z = 0.0 m. The numerical grid is oriented such that the gravity force is directed downwards and the yaxis is horizontal (with respect to the gravity). At the same time, the x-axis forms an angle of 15 degrees with respect to the horizon, and indeed, the gravity force forms an angle of 15 degrees with respect to the z-axis. Consequently, that originates a pressure gradient, and the flow goes toward the negative direction of the x-axis (see the blue arrows directed to the lefthand side).

Migration and distribution of a continuous leak of NAPL
All  Table 1 gives the material properties and parameter details used in the numerical simulations performed with CactusHydro [26] all over this work. In particular, it shows the density of the contaminant, which is 881kg/m 3 for the LNAPL, and 1200kg/m 3 for the DNAPL (see next subsection). It is important to stress that the numerical simulations with LNAPL and DNAPL differ only by their density values. The assigned input parameters used in the numerical simulations model have been taken from the Literature (not by conducting an experiment). They have been chosen to be representative of an aquifer in a "Loamy Sand" geological structure. The value of the density of 881kg/m 3 for the LNAPL corresponds to a value compatible with Crude Oil, while 1200kg/m 3 for the DNAPL was chosen in such a way to represent a generic, not too dense, one.
The porous medium is fixed to a value of porosity of 0.3. The absolute permeability is assumed to be k = 4.14 × 10 −10 m 2 all over the grid, with k x = k y = k z . The relative permeabilities and capillary pressure were obtained using Eqs (13)- (18). At zero saturation, the capillary pressure air-water is 793.80 Pa or, equivalently, 0.081 m, while the capillary pressure air-nonaqueous at zero saturation is 556.66 Pa or, equivalently, 0.0566 m.
After being released into the environment, the LNAPL starts to migrate downward the unsaturated zone under the influence of gravity (see Fig 2, the left-hand side, where the saturation (σ n ) contours are shown for different times). It moves predominantly downward. Lateral spreading may also occur due to the effect of the capillary pressures. Fig 2 second row, the lefthand side, shows the impact of the contaminant after 5.7 hours. The contaminant arrives at the groundwater table after two days and 8.9 hours (third row, the left-hand side), where a capillary pressure between air-contaminant and contaminant-water is present. Since the contaminant has a density lighter than the water (ρ n = 881kg/m 3 , see Table 1), it remains floating in the groundwater table zone representing a physical barrier. Part of the contaminant remains in the capillary fringe, part of it moves to the left direction due to a pressure gradient created by the gravity force's oblique position (see Fig 2 fourth row, left-hand side). Fig 2 also shows a few streamlines (in blue). Initially, they move undisturbed toward the left direction (left-hand side of Fig 2). They change their direction slightly under the contaminant's presence, which initially goes a bit down the groundwater table since the spill is sufficiently large, being a continuous leak source. On the right-hand side of Fig 2, the saturation contours are viewed in the y − x plane. The transient numerical simulation shows the behavior of this LNAPL after 9 days and 11.6 hours, although it is possible to go further in the numerical simulation. The first two rows (on the right-hand side) show no contaminant since the plane is located at z = 0 m (groundwater table). Only when the contaminant reaches this level may it be noticed. See the third row on the right-hand side how the contaminant arrives at the groundwater table, keeping its initial square-like section. After a while, this becomes an oval-like shape (third and fourth row, righthand side).  saturation is S n = 0.9 (red points). Below z = 8.0 m, this saturation goes abruptly to zero since initially the contaminant is situated on top of the parallelepiped (see also Fig 2, top, left-hand side). The water saturation (blue points), instead, is zero in the unsaturated zone (being completely dry) and is equal to one from z = 0.0 m to bellow (saturated zone). The air saturation (green points) equals 0.1 in the unsaturated region where there is also the contaminant, increases to one between z < 8.0 m and the groundwater table, and goes to zero in the saturated zone (below z = 0.0 m). Notice that the sum of the three phases' saturation is always one at any depth value. After two days and 8.9 hours, the LNAPL continuous leak source arrives at the groundwater table (see Fig 2 third row, left-hand side). See Fig 3 (middle), where the saturation of the LNAPL is kept constant to 0.9 (being a constant leak). It then arrives at the groundwater table (see red points), where it increases to 1.0 since the contaminant accumulates around z = 0.0 m before starting to move to the left-hand side together with the flow, to finally decrease very rapidly after passing through to the saturated zone where the sharp front goes to zero. Notice that the air saturation is now 0.1 above the groundwater table and zero in the saturated zone. After nine days and 11.6 hours (Fig 3 right-hand side, together with Fig 2 fourth row, left-hand side), the LNAPL remains partly in the capillary fringe with saturation one, partly goes into the saturated zone, being a continuous leak source (red points). Fig 4 shows three-dimensional numerical results on the pressure for the LNAPL (continuous leak source) as a function of the depth at different times. Initially, at time equal to zero (blue points), the pressure is equal to the atmospheric one (101325Pa), in the unsaturated zone composed by air-contaminant. Then the pressure increases, when moving from the groundwater table to the bottom, up to a value of 220 KPa. Similar behavior can be observed after two days and 8.9 hours and nine days and 11.6 hours (see the orange and green points, respectively). It can be noticed that the green points correspond to the contaminant that has already reached the groundwater table (see Fig 2). The pressure also increases slightly in the unsaturated zone.
3.1.2 DNAPL numerical results. Fig 5 shows similar numerical simulation results using the same set of values as the previous one (with LNAPL), except that the contaminant is a DNAPL with a density of 1200kg/m 3 (see Table 1 for details). The right-hand side shows the saturation contours in the (z − y) one. As in the previous case, the DNAPL is constantly released with a spill rate of 2.14 kg/s in the unsaturated zone, being a continuous leak source. It goes downward due to gravity, first in the unsaturated dry zone (see Fig 5, first and second row, that show this movement at time zero, and after one day and 4.4 hours). Notice how the fluid flows in the left direction (Fig 5 left-hand side) due to a pressure gradient. Instead, the right-hand side shows the z-y plane with a zero-gravity component effect and, thus, no privileged direction in the y − axis. That is the reason why the righthand side of Fig 5 is symmetric around the y − axis. When the contaminant arrives at the groundwater table, which acts like a physical barrier created by the different phases, unlike the previous case with an LNAPL, it keeps moving to the bottom (aquiclude) of the saturated zone (due to its density) while moving in the left direction due to the pressure gradient. Notice how the DNAPL arrives at −10 m depth after four days and a few hours and advances more rapidly on the left-hand side to the previous LNAPL case (see Fig 5 the third and fourth row, left-hand side). Fig 6 shows the three-dimensional numerical simulation results on the depth as a function of the water saturation S w (blue points), DNAPL saturation S n (red points), and air saturation S a (green points) at various times, for a continuous leak of DNAPL of Fig 5. The situation at zero time shows essentially the same behavior as the previous example with LNAPL (Fig 3). The contaminant saturation (red points) equals 0.9 for z greater than 8.0 m while it is zero in the rest of the unsaturated/saturated zone. The water saturation is different from zero only in the saturated zone, and the air saturation is 0.1 for z > 8.0 m, 1 for 0 m < z < 8.0 m, and zero elsewhere. After being released, the contaminant starts to migrate downward. This is represented in the middlebox of Fig 6, where the red points indicate the sharp front of DNAPL saturation moving along the vertical direction. The difference between this figure and the previous example with LNAPL (Fig 3) is that when the contaminant arrives at the groundwater table (see right-hand side of Fig 6), it keeps moving deeper into the saturated zone (that is why the red points are more immersed into the saturated zone). Notice how the sum of the three-phase saturations is always one (for a fixed depth value) and how the sharp contaminant front is immersed into the aquifer. Fig 7 shows three-dimensional numerical simulation results on the pressure as a function of the depth at different times, which is essentially similar to the previous case (in Fig 4), except that the green points here deviate from the other two (orange and blue) deeper than those in Fig 4.

DNAPL numerical results.
Consider now a small volume of leak contaminant in the unsaturated zone, a DNAPL, as is shown in Fig 8. The difference from the previous case is that the contaminant is now a finite volume (the other parameters are the same, see Table 1).    symmetric around the y − axis (right-hand side). Notice how the top of the contaminant begins to empty (the DNAPL saturation contours start to change from red to blue). When it arrives at the groundwater table, it keeps going to the bottom of the saturated zone since the DNAPL is denser than the water. See the third row of Fig 8. Contextually, the top of the initial position of the contaminant is almost empty, as can be seen from its saturation contour value. The difference between this situation and that of Fig 5 is that,  Notice how the air saturation is one where there is no contaminant. Later, the contaminant moves downward, and the DNAPL saturation decreases (center). Contextually, the air phase begins to occupy the space left by the contaminant (green points). Finally, (right-hand side), the contaminant arrives at the groundwater table and enters into the saturated zone replacing the water (see the last row of Fig 8, same time). Also, it can be observed (red points) how the contaminant in the upper zone is emptied. It can be clearly seen the shock front and the rarefaction of the DNAPL saturation in the unsaturated dry zone.

LNAPL numerical results.
Consider now a small volume of contaminant, an LNAPL, with a density of 881kg/m 3 (see Table 1 for details). Compared to the previous case, only the nonaqueous phase density changes.  The same situation is corroborated in Fig 11, which shows numerical simulation results of the depth as a function of the water saturation S w (blue points), LNAPL saturation S n (red points), and air saturation S a (green points) at various times for a small leak of LNAPL of Fig  10. Substantially, the principal difference between Figs 9 and 11 is that for the DNAPL, the contaminant (red points) goes through the saturated zone (bellow the groundwater table) for the LNAPL does not happen; therefore the saturation empties less.

Migration and distribution of a continuous leak of DNAPL in presence of impermeable "lenses" in heterogeneous aquifers
3.3.1 Impermeable "lens" in the unsaturated zone. Consider now the case in Fig 12 where a continuous leak source of DNAPL (similar to the case represented in Fig 5) encounters an impermeable zone (such as a clay lens in heterogeneous aquifer media), here simplified as a parallelepiped with an absolute permeability of 4.14 × 10 −14 m 2 (ten thousand smaller than the value of the rest of the domain, see Table 1). The (green) impermeable parallelepiped is sit-  Fig 5). After two days and 8.9 hours, the contaminant moves downward, and its saturation (red points) increases from 0.9 up to one, just above the impermeable at z = 6.0 m, as it accumulates in that area (center figure). Then it saturation abruptly goes to zero since we are showing the plane x = 0, and instead, the contaminant moves on the left-hand side (following the groundwater flow). This situation remains almost invariable for later times (right-hand side). x and z − y are shown at different times. Since the impermeable parallelepiped is a smaller one, with respect to the previous case, the continuous source of DNAPL has the opportunity to reach the groundwater table zone faster than the previous case and also from the other side of the parallelepiped situated at x = 5 m (see Fig 14, second row). Once arrived at the groundwater table/capillary fringe, the DNAPL keeps going downward and moving to the left direction due to a pressure gradient, as in the previous cases (see the third row). The last row of Fig 14 shows the saturation contour of the contaminant after seven days and 2.7 hours. Both concentrations situated below the parallelepiped are moving together toward the left side. Fig 15 shows the three-dimensional numerical simulation results of a depth as a function of the water saturation S w (blue points), DNAPL saturation S n (red points), and air saturation S a (green points) at various times for a continuous leak of DNAPL of Fig 14 (plane x = 0). The continuous source of DNAPL encounters an impermeable lens (depicted in green), as shown The rest of the unsaturated zone is filled with air (green points), while the saturated zone is filled with water (blue points). At later times, three days and 13.3 hours, the contaminant has already arrived at the impermeable lens and accumulates on top of it (middle). That is why the saturation sharp front increases from 0.9 to 1.0 (red points) and then goes abruptly to zero (in the region where there is the impermeable). On the right-hand side, there is the DNAPL saturation (red points) different from zero just below the groundwater table and is the DNAPL that passes under the impermeable in the plane x = 0. See Fig 14 fourth-row.  Fig 14). The figure was generated using an open-source post-processing named VisIt (https://wci.llnl.gov/simulation/computer-codes/visit). Notice how the contaminant moves on the impermeable parallelepiped zone and spills around the impermeable zone. The figure shows the iso-surface of equal contaminant density distribution after 8.15 days of the evolution shown in Fig 14. 3.3.2 Impermeable "lens" in the saturated zone. Suppose that the impermeable obstacle is situated in the saturated zone (instead of the unsaturated one). See Fig 17 (see the second  row), the DNAPL keeps going downward and moving to the left direction, due to a pressure gradient (see the third row), and eventually will reach the bottom zone being denser than the water (fourth row). Three-dimensional numerical simulation results of a depth as a function of the water saturation S w (blue points), DNAPL saturation S n (red points), and air saturation S a (green points) at various times for a continuous leak of DNAPL of Fig 12 (plane x = 0). The DNAPL leak encounters an impermeable obstacle (depicted in green in Fig 12). The left-hand side shows a sharp front of contaminant saturation (red points) situated in the unsaturated zone at z = [8.0,10.0] m. At later times, (center) after two days and 8.9 hours, the contaminant saturation increases to 1.0 due to an accumulation on top of the impermeable zone while part of it reaches the saturated zone. This situation remains almost invariable for later times (right-hand side) since we plot the x = 0 plane. https://doi.org/10.1371/journal.pone.0266486.g013 Fig 18 shows the numerical simulation results of a depth as a function of the water saturation S w (blue points), DNAPL saturation S n (red points), and air saturation S a (green points) at various times for a continuous leak source of DNAPL of Fig 17 (plane x = 0). The DNAPL encounters an impermeable lens (depicted in green) in the saturated zone, as shown in Fig 17. The left-hand side shows the situation after 5.7 hours where the continuous source of DNAPL is situated on top of the parallelepiped, similar to Fig 15. After two days and 8.9 hours (middle), the contaminant arrives at the groundwater table (see Fig 17, second row). Its saturation increases from 0.9 to 1.0 since it arrives at the saturated zone and tends to accumulate in that region, but at the same time keeps moving to the bottom (see red points), being denser than the water. Since the impermeable lens is situated just below the groundwater table, the contaminant saturation first decreases and then increases again. It finally goes to zero. This situation is better represented on the right-hand side after seven days and 2.7 hours. Three-dimensional numerical simulation results of a depth as a function of the water saturation S w (blue points), DNAPL saturation S n (red points), and air saturation S a (green points) at various times for a continuous leak of DNAPL of Fig 14 (plane x = 0). The DNAPL leak encounters an impermeable obstacle (depicted in green in Fig 14). After 5.7 hours, there is a sharp front of contaminant saturation (red points) situated in the unsaturated zone at z = [8.0,10.0] m. The DNAPL saturation increases from 0.9 to one at later times since it accumulates on top of the impermeable parallelepiped. Then abruptly goes to zero and finally increases in the saturated zone (righthand side). Notice that the DNAPL and water saturation sum is one for a fixed depth value. https://doi.org/10.1371/journal.pone.0266486.g015

Conclusions
In this work, we presented high-resolution three-dimensional numerical modeling investigation of the migration of three-phase immiscible fluid flow in variably saturated zones. We investigate the temporal evolution of the migration of the immiscible contaminant problem in a porous medium following the saturation contour profiles of the three-phases fluids flow. We considered both light nonaqueous phase liquid and dense nonaqueous phase liquid, initially released in unsaturated dry soil, and investigated several initial conditions, including impermeable parallelepipeds that mimic clay "lenses" in heterogeneous aquifer systems. The numerical simulations were obtained using CactusHydro code [18], based on the Cactus toolkit [36,37], that uses a high-resolution shock-capturing conservative method that precisely follows the advective part of the fluid flow. The results were validated through a classical convergence test running the same code at different resolutions. We show that it is possible to follow with high precision the migration of a contaminant sharp front (LNAPL or DNAPL, continuous or small volume source) in a variably saturated zone as a whole using a unique mathematical model which includes the complete set of mathematical equations expressed in terms of the saturation, permeability, capillary pressure, density, the viscosity of the three-phases. We show that the difference between the fate of a DNAPL and LNAPL (when the other parameters are the same) is just the density of the contaminant. We also show the numerical results of the three-phase saturation very precisely as a function of the depth. Next step of this research will be to compare numerical simulations results with laboratory experimental results such as a contaminant leakage in a sand tank.